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We illustrate how the systematic inclusion of multi-spin correlations of the quantum spin-lattice 
systems can be efficiently implemented within the framework of the coupled-cluster method by 
C"| examining the ground-state properties of both the square-lattice and the frustrated triangular- 

es ' lattice quantum antiferromagnets. Various physical quantities including the ground-state energy, 

the anisotropy susceptibility, and the sublattice magnetisation are calculated and compared with 
those obtained from such other methods as series expansions and quantum Monte Carlo simulations. 



^ ■ I. INTRODUCTION 

The techniques now available in the field of ab initio quantum many-body theory have become increasingly refined 
over the last decade or so. This is particularly true for whati-js-inowadays recognised as one of the most powerful 
1 modern techniques, namely, the coupled-cluster method (CCM)liTB. The results obtained from the CCM have become 
fully competitive with series expansions, variational calculations and quantum Monte Carlo (QMC) simulations (for 
q ■ the cases in which QMC may be applied). 

Quantum magnets not only provide useful models of many physically realisable magnetic systems but also serve 
as prototypical models of quantum many-body systems. Their rich phase diagrams due to strong quantum effects 
have naturally provided an excellent test-bed where the above-mentioned methods can be applied and furthtp-refined. 
, One example demonstrating rich and initially unexpected behaviour is provided by the Haldane conjectureLL3, which 
CNj ' states that the one-dimensional (ID) spin-1 Heisenberg antiferromagnet (HAF) possesses an excitation gap, in sharp 
contrast to its spin-i counterpart. This was surprising at the time because conventional spin- wave theory predicts 
a gapless excitation spectrum regardless of spin magnitude. However, the Haldane conjecture has subsequently been 
confirmed by numerical calculationsEM Moreover, in the aftermath of the discovery of the superconducting cuprates, 
much effort has been devoted to uncovering such subtle effects as spin-nematic, spin-Peierls and chiral spin liquid 
orderings in two-dimensional (2D) quantum antiferromagnets, among which the frustrated quantum aatiferromagnets 
on the triangular and the Kagome lattices have recently attracted considerable theoretical attentiorJl3i2l. 

The CCM has been applied to various quantum magnets-over the past six years. The first application of the CCM to 
these systems was performed by Roger and HetheringtonEa, who obtained good results at low levels of approximation 
for the ground-state energy of both the ID chain and the 2D square-lattice HAF, and also for solid 3 He where ring 
exchanges of nuclear spins are considered. Since then the CCM has been applied to thej-isotropic (Heisenberg) and 
anisotropic HAF (or XXZ model) in ID and on the 2D snuare lattice, both for spin-iE3 c3 systems and higher spins 
\ systemsE3E3; to the spin-LHeisenberg-biquadratic modeled; and to sncb frustrated spin models as the spin-i J1-J2 (or 
Majumdar-Ghosh\-mode]EEl~E3 and the 2D triangttl&p, lattice HAFlj£-X ft has also been applied to the spin-1 easy- 
. £h ! plane ferromagnetic. Among these, Bishop et a/B3'E3 not only put forth several systematic localised approximation 
schemes to perform higher-order calculations yielding good results on the ground-state sublattice magnetisations and 
approximate excitation spectra, but also used an infinite-order, two-body approximation scheme to obtain evidence of 
a zero-temperature quantum phase transitions. Moreover, the systematic inclusion-jof spin-spin correlations based on 
a dimerised state has also been made possible within the framework of the CCMEIrEil, to study spin-Peierls ordering. 
This may provide a possibkvineoad to probe more subtle topological order in the absence of solid order, as in the 
case of the chiral spin liquidli3 : tll. The quantitative description of such phases remains one of the most challenging 
problems for modern microscopic quantum many-body theory in general, and the CCM in particular. 

More recently, attention has been given to extending the CCM calculations to higher orderscJ in the particular 
case of the XXZ model, by using a localised approximation scheme, and by taking into account multi-spin correlations 
on up to 10 contiguous lattice sites in ID and on up to 6 contiguous lattice sites in 2D. The ground-state energies, 
for example, are found to be in excellent agreement (i.e., within about 0.03%) with the exact result in ID, and with 
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those obtained from spin-wave theoryEZI, series expansio rS and QMC calculations^ in 2D. However, it is fair to 
say that in 2D, to achieve the same accuracy on other more interesting physical quantities such as the ground-state 
sublattice magnetisation and the excitation spectrum, and to further clarify the nature of zero-temperature quantum 
phase transitions, the inclusion of multi-spin correlations of still higher orders is clearly needed. Since the extent of 
the task of determining the CCM equations and solving them grows extremely rapidly with the approximation level, 
efficient algorithms for performing the CCM calculations thus become indispensablea. 

The motivation of the present work is two-fold: (1) we analyse in detail the computational aspects of the CCM in the 
context of quantum antiferromagnets, in order to devise an efficient implementation and to show how the systematic 
inclusion of multi-spin correlations can be made simple; and (2) we revisit the spin-i quantum antiferromagnets on 
both the square and the triangular lattices by utilising this algorithm in a way which now enables us to increase the 
number of independent and localised multi-spin configurations to be considered by at least an order of magnitude 
over previous calculations. In this article we focus on the above two models in the regimes where a Neel-like order 
represents the corresponding classical limit. The present method, however, should be of general utility to quantum 
magnets where a spin-Peierls order is relevant, for example. 

A brief description of the contents of this article now follows. In Sec. II we describe how a reformulation of the CCM 
problem has been achieved. The characteristic CCM similarity transform of the spin operators is evaluated at a very 
general approximation level, and the Hamiltonian is reformulated purely in terms of spin-raising operators for a spin-i 
system. We also show that the form of the new Hamiltonian easily lends itself to a localised set of approximations 
called the LSUBm scheme. In so doing the CCM technique is itself clarified. The computational method used to 
determine our fundamental set of configurations within this approximation scheme is described, and the derivation 
of the resulting CCM equations is discussed in detail. The method is then anplied to the square-lattice spin-i XXZ 
antiferromagnet, and the triangular-lattice spin-j anisotropic antiferromagnettJ in Sees. Ill and IV respectively. Our 
conclusions are given in Sec. V, where we also discuss possibilities for further extending our computational solution 
to even higher-order calculations by making use of parallel processing and other strategies. 



II. THE CCM FORMALISM FOR SPIN-LATTICE MODELS 



A. Basic Ingredients of the CCM 

Since detailed descriptions of the fundamentals of the CCM are available in the literatureEHl, we only highlight the 
essential ingredients of its application here. The exact ket and bra ground-state energy eigenvectors, \*f?) and (^1, of 
a many-body system described by a Hamiltonian H , 

H\*)=E g \*); {9\H = E B {9\, (1) 

are parametrised within the single-reference CCM as follows: 

|*) = e 5 |<I>}; S = J2 SI C+, 

(*| = (<S>\Se~ s ; S = l + Y,hCj. (2) 

The single model or reference state |<E>) is required to have the property of being a cyclic vector with respect to 
two well-defined Abelian subalgebras of multi-configurational creation operators {C/} and their Hermitian-adjoint 
destruction counterparts {CJ = (C/)^}. Thus, |<I>) plays the role of a vacuum state with respect to a suitable set of 
(mutually commuting) many-body creation operators {Cf}, 

c/l*) = ,1^0, (3) 

with Cq = 1, the identity operator. These operators are complete in the many-body Hilbert (or Fock) space, 

i = |$)<$| + ^c+|$)($|C7. (4) 

1*0 

Also, the correlation operator S is decomposed entirely in terms of these creation operators {Cj~}, which, when 
acting on the model state ({C/|<I>)}), create excitations about the model state. We note that although the manifest 
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Hermiticity, ((*|* = |*)/(*|*)), is lost, the intermediate normalisation condition (^l*) = ($1*) = ($1$) = 1 is 
explicitly imposed. The correlation coefficients {sj,s/} are regarded as being independent variables, even though 
formally we have the relation, 

W S = ^ 5t si^ ■ ( 5 ) 
(<J>|e ST e s |<I>) 

The full set {s/,Sj} thus provides a complete description of the ground state. For instance, an arbitrary operator A 
will have a ground-state expectation value given as, 

A ee (*|A|*> = ($|Se- s Ae s |$) = A ({sj, §/}) . (6) 

We note that the exponentiated form of the ground-state CCM parametrisation of Eq. (||) ensures the correct 
counting of the independent and excited correlated many-body clusters with respect to |$) which are present in the 
exact ground state |\&). It also ensures the exact incorporation of the Goldstone linked-cluster theorem, which itself 
guarantees the size-extensivity of all relevant extensive physical .quantities. One crucial difference between jthe CCM 
parametrisation of the ground state and those used in spin-wavetiZl and variational Monte Carlo calculations^ is that 
although they all adopt an exponentiated form, the former (CCM) contains spin-raising operators only. 

The determination of the correlation coefficients {s/,s/} is achieved by taking appropriate projections onto the 
ground-state Schrodinger equations of Eq. (|l|). Equivalently, they may be determined variationally by requiring 
the ground-state energy expectation functional H({si,§i}), defined as in Eq. (^), to be stationary with respect to 
variations in each of the (independent) variables of the full set. We thereby easily derive the following coupled set of 
equations, 

8H/5~ Sl = 0^ ($|C7e- s iJe s |$) = 0, J^O; (7) 
6H/Ssi = 0^- ($\Sc- s [H,C+}e s \<S>) = 0, 1^0. (8) 

Equation (Q) also shows that the ground-state energy at the stationary point has the simple form 

E g = E g ({sj}) = <$ \c- s He s \$) . (9) 

It is important to realize that this (bi-)variational formulation does not lead to an upper bound for E g when the 
summations for S and S in Eq. (|^) are truncated, due to the lack of exact Hermiticity when such approximations are 
made. However, it is clear that the important Hellmann-Feynman theorem is preserved in all such approximations. 

We also note that Eq. (0) represents a coupled set of nonlinear polynomial equations for the c-number correlation 
coefficients {s/}. The nested commutator expansion of the similarity-transformed Hamiltonian, 

H = e- s He s = H + [H,S] + ^[[H,S\,S\ + ■ ■ ■ , (10) 

together with the fact that all of the individual components of 5* in the sum in Eq. (H) commute with one another, 
imply that each element of S in Eq. (||) is linked directly to the Hamiltonian in each of the terms in Eq. ([!(]). Thus, 
each of the coupled equations ([?]) is of linked cluster type. Furthermore, each of these equations is of finite length when 
expanded, since the otherwise infinite series of Eq. ( [To| ) will always terminate at a finite order, provided (as is usually 
the case) that each term in the second- quantised form of the Hamiltonian H contains a finite number of single-body 
destruction operators, defined with respect to the reference (vacuum) state |$). Therefore, the CCM parametrisation 
naturally leads to a workable scheme which can be efficiently implemented computationally. It is also important to 
note that at the heart of the CCM lies a similarity transformation, in contrast with the unitary transformation in a 
standard variational formulation in which the bra state (^| is simply taken as the explicit Hcrmitian adjoint of 



B. Computational Aspects of the CCM for Spin-Lattice Models 



1. Ket-State CCM Equations 

In this section, the general formalism of the CCM outlined in the previous section will be henceforth applied to 
spin-i quantum magnets with the emphasis on the common computational aspects involved in its implementation. 
Further model-specific details are deferred until Sees. Ill and IV. 
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To be specific, in this article we restrict ourselves to spin-^ quantum antiferromagnets in the regimes where the 
corresponding classical limit is described by a generalised Neel-like ordering, i.e., where all spins on each sublattice 
are separately aligned in the coordinates of a global quantisation axis. However, it is a simple task (see Sees. Ill 
and IV for details) to introduce a different local quantisation axis on each sublattice by a suitable spin-rotation 
transformation, such that the above Neel-like state becomes a fully aligned ("ferromagnetic") configuration in the 
local spin coordinates. This "ferromagnetic" state, |$), will be chosen as the uncorrelated CCM model state, where, 
in the local axes, all spins point along the respective negative z-axis, 

N 

|$) = | l)i ; in the local quantization axes . (11) 
i=i 

The correlation operator S is then decomposed wholly in terms of sums of products of single spin-raising operators, 
s% = sf + is v kl again defined with respect to the local quantisation axes, 

S=[i 1 }sf + [i 1 i 2 ]sfsf+--- , (12) 

where [ii], [1112] and so on stand for the corresponding (symmetric) spin-correlation coefficients (recall {s/} in Sec. 
IIA) specified by the sets of site indices, {«i}, 12} and so on, on the regular lattices under consideration. Implicit 
summations over repeated indices are also assumed. According to Eq. (^) , the spin-correlation coefficients in Eq. ( |T^ ) 
are to be determined by a set of CCM nonlinear equations: 

= (^K^---sJ M e- s He s \^) , (13) 
where s~ sj ■ ■ ■ sj is the Hermitian conjugate of the corresponding multi-spin correlation string sj sj~ • ■ • sj . 

31 32 3 m ■> ° fat- a 7, 7, 

In practice we clearly need an approximation scheme to truncate the expansion of S in Eq. ( 
infinite subset of the full set of multi-spin configurations {/}. The three most commonly used truncation methods up 
till now are: (1) the SUBn scheme, in which all correlations involving only n or fewer spins are retained, however far 
separated on the lattice; (2) the simpler SUBn-m sub-approximation, where only SUBn correlations spanning a range 
of no more than m adjacent lattice sites are retained; and (3) the systematic local LSUBm scheme, which includes 
all multi-spin correlations over all possible distinct locales on the lattice defined by m or fewer contiguous sites. Only 
the last approximation scheme will be adopted throughout this article. 

The first step in the practical implementation of the LSUBm CCM is to enumerate all of the distinct multi-spin 
configurations or correlated clusters, which we shall henceforth call fundamental configurations, 12, ■ ■ •■>in} with 
n < m, retained in the LSUBm approximation. It should be noted that the multi-spin configurations that are related 
by Hamiltonian symmetries, translational and rotational alike, are counted as one single distinct configuration. Such 
a correlated cluster can be either a connected cluster of size m (also called a "lattice animal" or "polyomino") or 
a subset of it (connected or disconnected) . Although the asymptotic behaviour of the number of lattice animals on 
a regular lattice remains an open combinatorial question, efficient algorithms for enumerating latticej-apimals up to 
sizes of about 20 have been developed in various fields including percolation and cell growth problemsc3. 

The second step in our modular implementation, namely, generating the corresponding set of CCM equations, 
is what we will focus on in the remainder of this section. Equation ( |l3| ) reveals that there are essentially two 
computational aspects involved in obtaining all possible non-zero contributions to its right-hand side. The first is 
to calculate the similarity-transformed Hamiltonian which then acts on the model state, and the second is to select 
terms of the similarity-transformed Hamiltonian that match exactly the string of spin-lowering operators represented 
by the set of site indices {ji, J2, ■ ■ ■ , Jm}- The first aspect is intrinsically related to the noncommutative nature of 
quantum spin operators, and the second to the geometric considerations of the lattice on which the Hamiltonian is 
defined. We address these two aspects in more detail below. 

The computation of the similarity-transformed Hamiltonian, H = c~ s He s , which acts on the model state |$) can 
be performed straightforwardly by making use of the relations = and s z |$) = — The goal here is to 

completely eliminate s z and s~, and thus retain the creation operators only, by utilising the commutation relations 
of the spin operators, namely, [s^s 1 * 1 ] = is^ and [s _ ,s + ] = — 2s z . This greatly simplifies the matching problem in 
generating the CCM equations as discussed below. To this end, we note that the similarity-transformed single-spin 
operators can be expressed as: 

4 = 4 , % = 4 + F *4 , $k = - 2F„4 - {F k fs+ , (14) 

where Fk = 53/ ■ ■ ■ i;_i]s^ ■ • • s+ . Furthermore, the commutation relations between the spin operators and the 
Fk operators can also be written in the following compact forms, 



12) to some finite or 
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[s z k , (F m ) 2 ] = 2F m G km s+ , [s k ,(F m ) 2 ] = -2(G km ) 2 s+-4F m G km s z k , (15) 

where G km = ~ l)[kmii ■ ■ ■ «/_2]s^ • • • sj_ 2 - Unlike in previous equations, repeated indices in Eq. ( |l5| ) do 

not imply summations. Clearly both F and G operators contain creation operators only. Consider then a typical 
two-spin interaction term, s^s^, for example, as contained in most spin- lattice Hamiltonians (see Sec. Ill and IV for 
a full description of the quantum spin Hamiltonians actually studied in this article). It is easy to prove the following 
relation: 

s k s m \<S>) = (2(G km ) 2 s+s++mF m G kmS +s+ + (F k ) 2 (F m ) 2 s+s+)\<S>) 
- (2G km F m s+ + F k (F m ) 2 s+ + (F k ) 2 F m s+ + 2G fem F fc s+) |$) 

+ (G km +F k F m ^) . (16) 

InEq. d||), the resulting terms from s k s m are classified into three categories as explicitly containing both s k and 
s+, either s k or s^ nl and neither s k nor s m , respectively. The reason for such a classification will become clear when 
we consider the second aspect, namely, generating the CCM equations. Thus, the first case is the simplest of all 
three to deal with, since the site indices of all terms in the case, including both k and m, are completely fixed up to 
permutations by the target set {ji,j2, ■ ■ '>3m} according to Eq. (|lg|). Although, unlike in the first case, only one of 
k and m in the second case must lie within the set {ji,j2, • • •, Jm}, the search for k or m in the matching problem 
can be easily performed once m or k is fixed. This comes about since the two-spin interactions with which we mostly 
deal are usually short-ranged. Typical examples are the nearest-neighbour interactions where k and m are simply the 
nearest neighbours as in the two models considered in this article. By contrast, neither index k nor index m in the 
last case must belong to the set {ji,j2, ■ ■ •, ju}- Nonetheless, for the LSUBra approximation scheme used here, both 
k and m must lie within a finite set of indices for which {ji, j 2 , ■ ■ •, j'm} is a subset. 

To be more specific, let us consider the term in Eq. (|l6|), F k F m , which can be written explicitly as: 

F k F m =J252( h + VQ* + i)^ 1 • • • ^Hroni ' ' ' ' ' ' << ' ' ' < 2 (17) 

h h 

where summation over repeated indices is implicitly assumed. Therefore, generating the part of the CCM equations 
due to this particular term F k F m amounts to determining all possible non-zero contributions to • • • sJ M F k F m |$) 

according to Eq. (|l7|). This is achieved by partitioning the target set {ji, j2, ■ ■ ^Jm} into two subsets {i\ ■ ■ ■ 
and {n\ ■ ■ ■ n; 2 } with l\ + 1% = M, followed by a search for the appropriate k and m in a nearby region that includes 
{i\ ■ ■ ■ and {n\ ■ ■ ■ ni 2 } respectively, such that both correlation coefficients [kii ■ ■ ■ and [mn\ ■ -^x^A-Aie the 
(symmetry-related) fundamental configurations of a given LSUBm approximation. Unlike earlier workE3c§E£l where 
the maximum number of fundamental configurations is limited to 100 or so, the present approach based on partition 
completely eliminates the costly procedure implemented previously for avoiding double occupancies of spin-i objects, 
and thus reduces the CPU usage a great deal. This optimal implementation becomes possible because all of the 
terms (e.g., F k F m ) in the similarity-transformed Hamiltonian (H) and, more importantly, their explicit structures are 
completely specified by the seemingly tedious reformulation of H in terms of F k , F m and Gkm operators. The CCM 
ket-state equations so obtained can then be solved by the standard Newton-Raphson method. 



2. Bra-State CCM Equations 

According to Eq. (||) in Sec. IIA, it is necessary to obtain both the ket-state correlation coefficients {s/} and 
the bra-state correlation coefficients {s/} in order to compute a general ground-state physical quantity such as the 
sublattice magnetisation. The task of generating the bra-state CCM equations (see Eq. (jg)) turns out to be simple. 
Firstly, note that this set of coupled equations is linear in {s/}, as is evident in Eq. (^|); secondly, a simple equality, 
5 2 H/SSj5si — S 2 H /SsiSSi, demonstrates that the bra-state equations can be readily generated from the already 
obtained CCM ket-state equations by appropriate differentiations. 

Similarly, in the context of spin-i quantum antiferromagnets, the S operator is in general decomposed entirely in 
terms of annihilation operators which are again defined with respect to local quantisation axes: 

S=l + \h]sr + [hh]sr s r +... , (18) 
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where [ii], \i1i2} and so on denote the corresponding bra-state spin-correlation coefficients (recall {sj} in Sec. IIA) 
specified by the sets of site indices {ii}, {11,12} and so on, which is the analogue of Eq. (|1J). Moreover, the same 
LSUBm truncation scheme is also adopted. This means that there exists a one-to-one mapping between the ket-state 
and the bra-state fundamental correlation coefficients which, to ease the burden of notation, are now denoted as {x r } 
and {x r } respectively. Now let SH/Sx r = P r {x\, x 2 , • • • , A) = denote the r th ket-state CCM equation, which is given 
in terms of the ket-state correlation coefficients and A, which stands for other parameters included in the Hamiltonian 
such as anisotropy, for example. Consequently H may now be written in terms of P r (xi, x 2 , • ■ ■ , A) and the bra-state 
correlation coefficients, {x r } as: 

H = P (xi, x 2 , ■ ■ ■ , A) + }^x r P r (xi,x 2 , ■ ■ ■ , A) , (19) 

r=l 

where Pq{x\,x 2 , • • • , A) denotes the zeroth-order CCM term, i.e., the ground-state energy expression, and Np is the 
number of fundamental configurations retained for a given LSUBm approximation level. Therefore, the s th bra-state 
equation may now be rewritten in terms of x r and P r {x\, x 2 , • • • , A): 

dP a {x 1 ,x 2 ,---,X) ^ dP r {x 1 ,x 2l - ■ ■ ,X) 

^ +2^- dx~ s =° > s=l,2,...,N F . (20) 

r— 1 

These linear equations may now be solved by using a standard decomposition technique for linearly dependent equa- 
tions, such as the LU decomposition method, once the ket-state correlation coefficients {x r } are known. 

Equipped with the above CCM algorithms, which are readily implemented for spin-^ quantum magnets, we report 
our new findings on ground-state properties of both the square- and triangular-lattice spin-i quantum antiferromagnets 
in Sees. Ill and IV, respectively, in order to illustrate the relative simplicity of the present approach. Similar 
algorithms have also been successfully implemented to study the excitation spectra. This will be the subject of a 
future publicationca. 



III. SPIN-i XXZ ANTIFERROMAGNET ON THE 2D SQUARE LATTICE 

In this section we shall consider the spin-i XXZ model on the infinite square lattice. The XXZ Hamiltonian is 
given by, 



H = 



E 



<4 + s\a) + Asts* 



(21) 



where the sum on runs over all nearest neighbour pairs and counts each pair only once. The square lattice XXZ 
model has no exact solution, unlike its ID counterpart, although approximate analytical and numerical calculations 
have been performed. To put later CCM calculations in context, we note that the XXZ model has three regimes: 
an Ising-like phase characterised by non-zero Neel order; a planar-like phase in which the spins in the ground-state 
wavefunction are believed to lie in the xy plane; and a ferromagnetic phase. A Monte Carlo study of the 2D anisotropic 
Heisenberg model was performed by Barnes et alE3. They observed that the staggered magnetisation in the z-direction 
is non-zero for A > 1, but then appears to become zero below A = 1. They therefore conclude that-the critical point 
is probably very near to this point. In contrast to this Monte Carlo calculation, Kubo and Kishiu3 have used sum 
rules to investigate the ground state of this system. They state that the ground state possesses an off-diagonal long- 
range order (LRO) akin to that of the XF-like state at small anisotropy, 0.0 < A < 0.13 . Also, for A > 1.78 they 
observe that the system demonstrates non-zero Ising-like LRO. At A = — 1 there is a first-order phase transition to 
the ferromagnetic phase for this model. 

The isotropic Heisenberg point has been extensively studied using various approximate methods, and so shall be 
used as a test case for the CCM results discussed later in this section. RungecJ has performed the most accurate 
Monte Carlo simulation to date. This provides a value for the ground-state energy per spin of —0.66934(4), and a 
value for the subiattice magnetisation which is 61.5%±0.5% of the classical value. In comparison, linear spin-wave 
theory (LSWT)EII gives a value of —0.658 for the ground-state energy, and a value for the subiattice magnetisation 
which is 60.6% of the classical value. 
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A. The Model State 



We begin the CCM treatment of this spin system by choosing a suitable model state |$) (for a particular regime), 
such that all other possible spin configurations may be obtained by the application of linear combinations of products 
of spin-raising operators to this state. In the Ising-like regime, characterised by non-zero Neel order, a natural choice 
for the model state is the Neel state with spins lying along the z axis. (For clarity, this state will be referred to as 
the z-axis Neel model state throughout this article.) We note however that this model state is not the best choice 
for all values of A because the ground-state wavefunction of the XXZ model in the region — 1 < A < 1 is believed to 
contain only spins which lie in the xy plane. In this regime, we again use the classical Neel state, but this time with 
spins lying on the x axis. (This state will be referred to as the planar model state throughout this article.) Hence, we 
see that even for the same spin model and lattice a different choice of model state may be preferable, depending on 
the regime that we are investigating. 

We shall consider the Ising-like regime first, and, so that spins on either sublattice may be treated equivalently, we 
peform a rotation of the local axes of the up-pointing spins by 180° about the y axis. The transformation is described 
by, 

s x -s x , s y s v , s z -> ~s z . (22) 

The model state now appears mathematically to consist of purely down-pointing spins which is precisely given by 
Eq. (|ll|), and the Hamiltonian may be written in terms of these local axes as, 



H 



N 

E 



44 



2A& 



(23) 



For the planar model state, we again rotate the local axes of these spins on the separate sublattices such that all spins 
appear to be lie along the negative z direction which is again given by Eq. (O) . This is achieved by rotating the axes 
of the left-pointing spins (i.e., those pointing along the negative x directionjin the planar model state by 90° about 
the y axis, and by rotating the axes of the right-pointing spins (i.e., those pointing along the positive x direction) by 
270° about the y axis. (The positive z-axis is defined to point directly upwards, and the positive cc-axis is defined to 
point directly to the right.) Hence the transformation of the local axes of the left-pointing spins is described by, 



and the transformation of the local axes of the right-pointing spins is described by, 

s — > — s , s u — > s" , s — > s ; 
The transformed Hamiltonian for the planar model state is now given by, 

(A + l)( s + s + + srsj) + (A - l)(sfsj + srsj) - 



H 



1 N r 

-Y 



(24) 



(25) 



(26) 



In the remainder of this section, the power and flexibility of our new formalism is illustrated by focussing primarily 
on the planar model state applied to the XXZ model on the square lattice. Note however that equivalent calculations 
have been undertaken for the z-axis Neel model state, and that a general explanation is presented here only. A more 
detailed explanation of CCM calculations that deal with both the ground-state properties and the-axcitation spectrum 
using the new formalism for the XXZ model based on this model state will be presented in Ref.c2l. 



B. Fundamental Excitation Configurations: Lattice Animals 



As described in Sec. IIB, the first step of our modular solution is to obtain the set of fundamental configurations 
for a given approximation scheme by utilising appropriate lattice symmetries. Another constructive way to define 
the LSUBm scheme used here for the square lattice is to consider a right-angled bounding triangle containing m 
lattice points along the sides parallel to the axes (see Fig. [I] for a diagram of this construction where the bounding 
triangle for LSUB4 is shown). All possible fundamental configurations for the LSUBm approximation are then 
confined by this bounding triangle. This comes about because it is easy to show that all connected configurations 
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of size m (or the lattice animals of size m) are constrained to lie within or on this bounding triangle. Furthermore, 
as first shown by LunnonCj, the introduction of this bounding triangle greatly simplifies the recursive procedure of 
growing a connected cluster of given size. The disconnected configurations for LSUBm scheme are then constructed 
by successively considering all "subsets" of each member of the fundamental set of connected configurations, and all 
possible disconnected configurations are thereby generated. (The "subsets" here refer to all independent configurations 
which are formed by removing one or more spins from these connected configurations.) 

To be specific, for the square lattice, there are four rotational operations, (0°, 90°, 180°, 270°), and four reflections, 
(along the x and y axes, and along the lines x=y and x=—y), which preserve the symmetries of both the lattice and 
the Hamiltonian. Moreover, the Hamiltonian of Eq. (|2^), which is defined with respect to the planar model state, 
contains only even products of spin- flip operators and a single term containing two s z operators. Repeated application 
of this Hamiltonian to this model state yields the ground state (assuming that this model state is not orthogonal 
to it). Therefore the ground state will have an even numbers of spin flips with respect to this model state, and so 
we restrict the LSUBm approximations to contain even numbers of spin-raising operators in the ket-state correlation 
operator, S, only for this planar model state. As an example, in Fig. |l| we show all 10 fundamental configurations 
retained in the LSUB4 approximation when the planar model state is used in the CCM calculation. 

Further reduction in the number of fundamental configurations can be made when the z-axis Neel model state 
is used in the CCM calculations. This comes about because, although the total uniform magnetisation Sj,=Y2 i s l 
(where sf is defined with respect to a global quantisation axis and the sum on the index i runs over all lattice 
sites) is always a good quantum number independent of the model state used, only the z-axis Neel model state is an 
eigenstate of the total uniform magnetisation Sj.. In contrast, the planar model state is not an eigenstate of the total 
uniform magnetisation s^. Therefore, for the z-axis model state case one can explicitly conserve Sj, by restricting the 
fundamental configurations to those which produce no change in sf, with respect to the z-axis Neel model state. This 
restriction, for example, reduces the number of the fundamental configurations retained in the LSUB4 approximation 
to 7 if the z-axis Neel model state is employed in the CCM calculations, and see Fig. |l|. We tabulate the number of 
of fundamental configurations up to the LSUB8 level of approximation for both model states in Table |. Note that 
only the CCM calculations based on the z-axis Neel model state are actually carried out up to LSUB8 approximation 
in this article. 

C. Similarity-transformed Hamiltonian and CCM Ket-State Equations 

In order to solve the Schrodinger equation of Eq. ([!]), we shall specifically utilise the Hamiltonian of Eq. j2^), 
although a comparable analysis can also be performed for Eq. d23). The expression for the ket-state correlation 



operator of Eq. ( J12| ) is now used to write the similarity-transformed Hamiltonian, H, of Eq. (10) in terms of the 
operators F k and G km . Furthermore, we subdivide H into three categories as discussed in Sec. IIB to clarify the 
problem of finding the CCM equations, where H\<&) = e~ s He s \<f>) = (H 1 + H 2 + H 3 )\§), such that: 



kp k 



#i = o > 'A -i G km + F k F m ) + i(A - \){Fl + F%) \ s+s 
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~ \( A + !) E{ 2G L + ±G km F k F m + F^s+s+ ; (27) 
^ = i E { F ™ s ™ + F * s t + \0- - A )( F m4 + 

kp 

+ ~(A + 1) E j(2Gfem + F k F m ){F m s+ + F k s+ 

kp ** 

J ff 3 --^E{ 1 + ( A + 1 )( Gfc ™ + FfcF '™)} • 



+ ) \ ; (28) 

(29) 



Note that k runs over all lattice sites and that m is given by m = k + p. such that p covers all nearest neighbours to 
k. Hence we see from Eq. (p9|) that the ground-state energy of the XXZ model for the planar model state is given by, 



N 



2ari(A + l) + l 



(30) 
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where x\ = [k, k+p] represents all nearest- neighbour, two-body correlation coefficients in Eq. ( [l2| ) which are equivalent 
under the translational and rotational symmetries of the lattice, and z represents the lattice coordination number 
(i.e., the number of nearest neighbours to a given site). Note that the ground-state energy of Eq. (|3(]) is exact in the 
sense that both the exact expansion for S and any non-trivial approximation of it will always produce this expression. 



D. Results 



1. Ground-State Energy 

The ground-state energy for the planar model state is illustrated by Fig. |[ and we can see that these results 
appear to be in good agreement with Monte Carlo resultsa. The highest approximation that has been attempted 
for the planar model state is the LSUB6 approximation, which contains 131 fundamental configurations and gives 
a ground-state energy per spin of —0.66700 at the Heisenberg point. For the z-axis Neel model state, due to the 
reduced number of fundamental configurations, we were able to solve the LSUB8 approximation, which contains 1287 
fundamental configurations and gives an energy per spin of —0.66817 at A = 1. Hence, by utilising the new formalism 
we have increased the number of fundamental configurations used in a CCM-calculation for the z-axis Neel model state 
by over an order of magnitude compared to the previous best (i.e., LSUB6c3 which contains 75 configurations). Table 
U summarises the information regarding numbers of configurations and ground-state energies at A = 1. Note that the 
calculations for both model states at the Heisenberg point give exactly the same results at equivalent approximation 
levels, and so only one figure for the ground-state energy is quoted in Table Q. The reason for the equivalence is 
that the Hamiltonians of Eqs. ( |2^ ) and ( |2^ ) become identical at A=l. Also, all of the correlation coefficients for 
configurations at a given LSUBm level contained in the planar model state case but not contained in the z-axis Neel 
model state case become identically zero at A — 1. 

The results for the z-axis Neel model state are found to be less accurate than those using the planar model state 
in the region — 1 < A < 1. Conversely, for A > 1, the results based on the z-axis Neel model state become the more 
accurate of the two sets of calculations. This therefore vindicates our decision to use two separate model states in 
order to investigate the Ising- and planar-like phases of this model. 



2. Anisotropy Susceptibility 

Beyond certain values of the anisotropy parameter (called critical points) it is found that there is no physically 
reasonable solution to the LSUBm CCM equations for m > 4. This characteristic breakdown of the solution CCM 
equations has previously been related to a phase transition of the real systemEZl. We may also define a quantity called 
the anisotropy susceptibility, given by, 

■ ■ - gWg) (31) 

which diverges at the LSUBm critical points. Table [j] illustrates two sets of estimates for the critical points for 
the XXZ model based on the planar model state, corresponding to the ferromagnetic and antiferromagnetic phase 
transition points. Encouragingly, the critical points corresponding to the ferromagnetic phase transition become closer 
to A = — 1 with approximation level. Table [| also includes the estimates for the critical points obtained for the z-axis 
Neel model state corresponding to the antiferromagnetic phase transition point. Note that LSUBm results based on 
these two model states always bound the Heisenberg point, at which the true antiferromagnetic phase transition is 
believed to lie, and also appear to converge with increasing m. 



3. Sublattice Magnetisation 

We now consider a simple order parameter called the sublattice magnetisation, M + = — 2(s z ), which is defined in 
terms of the local, rotated spin axes. Hence M + is given by, 

_„ N No N F 

M+ = — £(*| s ?|*) = l__ ^\§F k s+^) = l-2j2(n r y.x r x r , (32) 

i=l k=l r=l 
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where the index i runs over all Nq sites of a sublattice, and Np again indicates the total number of configurations 
for a given LSUBm approximation level. Note that x r and x r , respectively, are the bra- and ket-state correlation 
coefficients associated with the r th fundamental configuration, and that n r is the number of spins in this configuration. 
From Eq. (|3^) we can see that in order to obtain a numerical value for the sublattice magnetisation we must first 
know the values of both the ket- and bra-state correlation coefficients. The manner in which we determine these 
coefficients was described in Sec. II. Our results for the XXZ model using the planar model state are shown in Fig. [|, 
and we note that these results appear to converge to a non-zero value as one increases the approximation level. Hence 
our results indicate non-zero Neel-type long-range order in the xy plane for the square lattice in the planar regime. 
Table | also summarises the results for the sublattice magnetisation at the Heisenberg point. Note that results based 
on both model states are again found to be identical at this point and so only a single number is quoted in Table |. 

In summary, the new CCM formalism has been used in order to calculate estimates of the ground-state energy 
and sublattice magnetisation for the 2D XXZ model, and these results are found to be in good agreement with 
other approximate calculations. The highest-order approximation for the z-axis Neel model state has been extended 
to LSUB8 level using the new formalism, which is an increase of over an order of magnitude in the number of 
fundamental configurations used in the previous-highest LSUB6 calculation. Also, the positions of the phase transition 
points obtained by the CCM, using both model states, are fully consistent with the known behaviour for this model. 
Finally, the results presented here also support the idea that this model contains both Ising- and planar-like phases. 



IV. SPIN-i TRIANGULAR LATTICE ANTIFERROMAGNET 



Unlike the square-lattice spin-^ HAF where various calculations including extensive quantum QMC simulation; 




strongly support the existence of a Neel ordering with a reduced magnetic moment of about 62% of its classical 
value (See Sec. Ill), the three-sublattice orderingJor the corresponding triangular case is much less clear. For 
instance, early variational wavefunction calculationsEil that include long-ranged two-spin and nearest-neighbour three- 
spin correlations support an ordered ground state with a value of sublattice magnetisation, M + = 0.68, i.e., as large 
as 68% of the classical value. Based on this antiferromagnetic correlated triaLyavefunction, fixed-node Green function 
Monte Carlo-simulations were recently performed on lattices of up to 324 sitesEI. These yielded a similar magnetisation, 
M+ fa 0.6(H However, series expansion calculationstJ, utilising up to llth-order terms in an Ising-like anisotropy 
parameter suggest that the triangular HAF may be at, or at least close to (with the magnetisation being extrapolated 
to a value of M + ps 0.20), the critical point of losing maepjetic order. This scenario has received support from exact 
diagonalisation calculations on lattices of up to 36 sitesc3. Yet another careful analysiso of the same data from 
exact diagonalisations in terms of a consistent description of the symmetries and dynamics of the quasi-degenerate 
joint states indicates the presence of sublattice magnetic-order with the magnetisation value, M + ss 0.50, which is 
also consistent with second-order spin-wave calculationsEll. Clearly, further work is still needed to account for the 
discrepancy, and to provide a more definite and converged result. 

Hence, in this section we therefore further apply the CCM to the spin-| triangular HAF and focus on model-specific 
details of the algorithm implementation discussed in Sec. II. Compared with earlier applications of the CCM on the 
triangular HAF which only include two-spin (though long-ranged) correlationsEij, the current calculations, which take 
into account all multi-spin correlations on up to six contiguous sites, have obtained various ground-state properties 
that are now found to be fully competitive with those obtained from the above-mentioned methods. 



A. Model Hamiltonian and Model State 



The spin~2 triangular HAF is described by the antiferromagnetic-coupling Hamiltonian, 

H = $i ■ > ( 33 ) 

where Sj denotes the spin-^ operator at site i on the infinite triangular lattice. The sum in Eq. (|33| ) on runs 
over all nearest-neighbour pairs and counts each pair once. We note that the operators in Eq. (p3|) are defined in 
terms of some global spin quantisation axes referring to all spins, whereas henceforth we shall consistently employ a 
notation in which the spin operators are described in terms of local (Neel-like) quantisation axes for each of the three 
sublattices (A, B, and C) of the triangular lattice. The classical ground-state of Eq. ( |33"| ) is the Neel-like state where 
all spins on each sublattice are separately aligned (all in the xz-plane, say). The spins on sublattice A are oriented 
along the negative z-axis, and spins on sublattices B and C are oriented at +120° and —120°, respectively, with respect 
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to the spins on sublattice A. In order both to facilitate the extension of the isotropic Heisenberg antiferromagnet to 
include an Ising-like anisotropy first introduced by Singh and Husea and to make a suitable choice of the CCM model 
state, we perform the following spin-rotation transformations. Specifically, we leave the spin axes on sublattice A 
unchanged, and we rotate about the y-axis the spin axes on sublattices B and C by —120° and +120° respectively, 



2 ° 2 ° ° 2 u 2 
1 z V3 x 



2 ^ " 2 S& ' (34) 

We may rewrite Eq. ( |33| ) in terms of spins defined in these local quantisation axes for the triangular lattice with a 
further introduction of an anisotropy parameter A for the non-Ising-like pieces, 



^(stsJ + srsD-^istsf + srsj)} , (35) 



H = 

+£k +s 7 + s * r 4)-| A 



where A = 1 corresponds to the isotropic Heisenberg Hamiltonian of Eq. (^) . We note that the summation in Eq. 
( |35| ) again runs over nearest-neighbour bonds, but now also with a directionality indicated by (i — > j), which goes from 
A to B, B to C, and C to A. When A = 0, the Hamiltonian in Eq. ( |35| ) describes the usual classical Ising system with 
a unique ground-state which is simply the fully aligned ( "ferromagnetic" ) configuration in the local spin coordinates 
described above. We choose this state as the uncorrelated CCM model state |$) which is, of course, precisely given 
by Eq. «. 



B. Fundamental Excitation Configurations: Lattice Animals 



Unlike the square lattice case discussed in Sec. IIIB where all lattice point-group symmetries are employed to 
produce symmetry-distinct configurations, care must be exercised here since not all of the lattice point-group symme- 
tries leave the lattice-spin Hamiltonian invariant. The Hamiltonian of Eq. (|3^) (or the CCM model state) explicitly 
breaks some of the lattice symmetries because of the presence of bond-directionality in the Hamiltonian. Thus only 
6 (instead of the full 12) point-group symmetries should be used in the symmetry reduction. These are, specifically, 
three rotational operations (0°, 120°, and 240°) together with three reflections about the lattice axes (i.e., lines that 
coincide with the edges of the triangular lattice). For example, the three configurations (a), (b), and (c) shown in 
Fig. |] are symmetry equivalent, as are the three configurations (d), (e), and (f). However, the former are regarded 
as inequivalent to the latter in the context of the present spin-lattice Hamiltonian problem. For the purpose of com- 
parison with the case in, say, percolation problems, and for the sake of concreteness, let us consider the connected 
configurations of size 6. If all 12 point-group symmetries were used, we would have obtained 82 symmetry-inequivalent 
configurations as shown in Fig. 5, which are further classified into two groups: 17 in group group A and 65 in group 
B. The configurations in group A are of higher symmetries than those in group B, and thus do not lead to new 
symmetry-inequivalent configurations when only 6 point-group symmetries must be used in the symmetry reduction 
as discussed above. Each configuration in group B, however, results in another new symmetry-inequivalent configura- 
tion. Therefore, the total number of symmetry-inequivalent connected configurations of size 6 is 147 for the present 
spin-Hamiltonian problem. The set of symmetry-inequivalent configurations (connected and disconnected) thereby 
forms the set of fundamental configurations. We tabulate the number of fundamental configurations up to the LSUB7 
level of approximation in Table ||. Note that only CCM computations up to the LSUB6 level of approximation are 
actually carried out in this article. 



C. Similarity- Transformed Hamiltonian and CCM Ket-State Equations 

Using Eq. ( |15[) given in Sec. IIB, we can straightforwardly carry out the similarity transformation of the Hamiltonian 
given by Eq. (pq); the resulting terms are further classified into three categories for reasons already discussed in Sec. 
IIB, i.e., H\$) = e- s He s \<P) = (H x + H 2 + H 3 )\$), as indicated below: 
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Hi = {- 2 ( G k>n + F k F m ) ~y + ^ X ( F * - + 2G fem + F k F m )\ s+s+ 

+ 7 E \-^ F "n + F ') - 3AGL ~ ^G km F k F m - Y F i F A "t'm ; 

kp ^ J 

+ \ 12 [^ X (G km + F k F m )(s+ - s+) + y (2G km + F k F m )(F m s+ + F k s+)^ 



^fe) — ^{G km + F k F m ) 



(36) 



(37) 
(38) 



Here the summation over k runs over all triangular lattice sites, while the summation over p is over the three directed 
nearest- neighbour vectors that point from A to B, B to C, and C to A as required by the explicit bond directionality 
in the Hamiltonian given by Eq. (|35|), and the index m = k + p. Each fundamental configuration for a given LSUBm 
approximation is then pattern-matched to each of the total 35 terms in the above similarity-transformed Hamiltonian 
following the procedure outlined in Sec. IIB to generate the entire set of coupled CCM ket-state equations. The CCM 
ket-state equations may then be solved by the Newton- Raphson method for nonlinear equations. Specifically, we start 
from the point A = 0, at which we know the exact solution where all the correlation coefficients are zero. We then 
use this known solution as an initial input in solving the CCM equations for a slightly increased nonzero anisotropy 
A. This procedure is carried out recursively to obtain the numerical results reported below. 



D. Results 



1. Ground-State Energy 



In Fig. |6| we show the ground-state energy per spin E g /N as a function of the anisotropy parameter A for various 
LSUBm approximations. The corresponding values at the isotropic Heisenberg point are also tabulated in Table 
0. The highest-order calculation, LSUB6, which consists of 758 independent fundamental correlation coefficients, 
yields E g /N — —0.54290. This value shpuld be compared with the value —0.5445 extrapolated from-finite-cluster 
diagonalisations of up to 36-spin clusterala, and the value — 0.5431±0.0001 from a recent QMC simulationcH. Compared 
with the corresponding classical value of —0.375, it is safe to say that the LSUB6 CCM calculation captures at least 
99% of the quantum corrections. Due to the lack of rigourous finite-size scaling results, the investigation of which 
itself presents an interesting and important open question, no attempt at extrapolation is made here. Nonetheless, 
the ground-state energy obtained from the LSUB6 approximation alone is fully competitive with those obtained from 
the above-mentioned alternative methods, and is clearly among the best estimates currently available. 

To make further contact with the highest-order series expansion kno wn t o dateEj, we have computed the perturbative 
solution of Eg/N in terms of the anisotropy parameter A. In Table III we tabulate the expansion coefficients from 
the LSUB6 approximation, together with the corresponding results from exact series expansionsEJ. We note that the 
LSUB6 approximation reproduces the exact series expansion up to the 6th order. This result lends further strong 
support to the conjecture that the LSUBm approximation reproduces the exact series expansion to the same mth 
ordered. Moreover, the fact that the corresponding values of several of the higher-order expansion coefficients from 
both the CCM LSUB6 perturbative solution and the exact series expansion remain close to each other shows that 
the exponential parametrisation of the CCM with the inclusion of multi-spin correlation up to certain order also 
captures the dominant contributions to correlations of a few higher orders in the series expansions. Further detailed 
series analysis, using such methods as Pade approximant techniques, for the perturbative-expansions of both the 
ground-state energy E g /N and the sublattice magnetisation M + will be reported elsewhereEll 



2. Anisotropy Susceptibility 



As already displayed in Fig. |6J, the LSUBm ground-state energy curve for m > 3 terminates at lower and upper 
critical values of the anisotropy, A Cl and A C2 respectively, beyond which no physical solution of the CCM ket-state 
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equations exists. We have also tabulated A C; and A C2 for various LSUBm approximations in Table Q. Although the 
CCM based on the model state given in Eq. ( Jll| ) with the three-sublattice magnetic ordering is bound to break down in 
the region of the anisotropy parameter space where the tuie,ground-state wavefunction possesses a different symmetry 
from that of the model state, it has been strongly arguedocJ that the terminating points may correspond to the critical 
points of a phase transition. In order to shed further light on the nature of the terminating points, we investigate 
the singular behaviour of the CCM correlation coefficients near these points by calculating the derivatives of these 
coefficients with respect to the anisotropy parameter. Specifically, we display in Fig. ^ the anisotropy susceptibility 
defined analogously to Eq. ( |3l|) as, Xa = —d 2 {E g /N)/d\ 2 . Clearly, Xa diverges at the corresponding terminating 
anisotropy parameters A Cl and A C2 . r-^ote that the lower terminating point A Cl clearly converges to a value of about 
—0.5, as argued by Singh and Huset£|. However, the upper terminating point A C2 , as tabulated in Table ||, remains 
considerably larger than the value of unity which was obtained via Pade analysis of the series expansion by Singh and 
Husaia. Further work is needed to explain this disagreement and also to determine the singularity exponents of the 
ground-state energy at the two critical points. 

3. Sublattice Magnetisation 

Once the ket- and bra-state correlation coefficients are known it is possible to evaluate the sublattice magnetisation, 
M + = — 2(s z ), which is similarly defined by Eq. ( |32| ) except that the subscript i now covers all N sites on the 
sublattice A, one of the three sublattices. In Table [il] we also tabulate the sublattice magnetisations at the Heisenberg 
point for the various LSUBm approximations. The highest-order LSUB6 approximation gives give rise ia a value of 
0.6561, which agrees well with the corresponding value of 0.60 obtained recently from QMC simulationso. However, 
the discrepancy with that from finite-size and the second-order spin- wave calculations still needs to be accounted for. 

The divergence in (s 2 ) seen in Fig. |^ near the critical points is a natural consequence of the approximate nature of the 
calculation. As we approach a critical point for a given LSUBm approximation, one of the CCM correlation coefficients 
x r becomes very large. The contribution to (s z ) from this coefficient also becomes very large and so the sublattice 
magnetisation diverges. The puzzling "upturn" of M + observed for the LSUB5 and the LSUB6 approximations near 
their respective upper critical points A C2 remains elusive to us at the present. 

In summary, compared with earlier applications of the CCM on the-triangular HAF which have already revealed 
interesting oscillatory behaviour in long-ranged two-spin correlational, the present high-order CCM calculations 
(with a systematic inclusion of multi-spin correlations on up to six contiguous lattice sites) have already obtained 
results including ground-state energy and sublattice magnetisation that are fully competitive with those obtained from 
other methods, and which are among the best available. Further detailed analysis of the large number of ket-state 
coefficients already obtained for the LSUB6 approximation, in order to achieve a better understanding of the nodal 
surface of the ground-state wavefunction, may provide a more microscopic justification of, and an extension to, the 
above-mentioned variational wavefunction. This may in turn lead to a better trial wavefunction for QMC simulations. 

V. CONCLUSIONS AND OUTLOOK 

In conclusion to this article, we restate the main points of the new formalism, and discuss the results for the two spin 
models to which it has been applied. We also consider the future development of the CCM using the new approach, 
and make suggestions regarding possible strategies in order to extend the CCM calculations to even higher orders. 

Our new formalism has reformulated the Hamiltonian purely in terms of spin-raising operators acting on some 
suitable model state. We note that the reduction in CPU processing time is realised by cutting out one large task, 
namely, the evaluation of the similarity transform. Also, the method presented here largely negates the task of checking 
for double occupancy for the spin-^ system, which can also cause considerable delay. The simple and straightforward 
nature of the terms within the similarity-transformed Hamiltonian means that one may easily amend the code in our 
programs to deal with other systems and lattices, and so we have a very flexible approach. 

We have seen that the ground-state energy results for both of the systems considered are fully competitive with 
QMC energies, without recourse to large computing facilities; we have used a HP series 700 workstation containing 
96 Megabytes of RAM to perform the calculations presented here. We have also seen that sublattice magnetisations 
have been determined which are again fully competitive with QMC calculations. 

The present limitations of the new method arise form two sources. One limitation is that we must save our CCM 
equations in memory or on disk in order to obtain expectation values. The rapid increase in the number of fundamental 
configurations with increasing LSUBm approximation level means that the cost in memory and disk space rapidly 
becomes prohibitive, and the task of solving the CCM equations becomes slower also. Another limitation is that, as 
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discussed in Sec. IIB, we search over a wide area in order to determine the 'free' indices within specific terms in H . 
In ID, we find that this area depends on m linearly, but in 2D, for example, this area increases quadratically with 
to. We need to perform this search because, as yet, we do not take into account the manner in which the members 
of the set of fundamental configurations are inter-related. For example, some configurations are related to others 
by the addition or subtraction of a single spin. Thus, most of the processing time is spent sweeping through this 
area regardless of whether we know that a particular configuration will contribute or not. A useful gain could be to 
determine the way in which the configurations are related, using graph theory^ for example, and then to use this 
knowledge in order to reduce the search area. 

We have seen that the CCM LSUB6 approximation fully reproduces the first 6 coefficients of the exact series 
expansion for the triangular HAF. Also, for higher orders in the series than sixth order, it is found that the CCM 
captures about 99% of the expansion coefficient of the 7th order, and about 75% of those of higher orders up to 11th 
order in the ground-state energy perturbation series on the triangular lattice. CCM expansions of the ground-state 
energy have also been performed up to sixtieth order, at a fraction of the computational cost compared to a direct 
calculation of the exact series to this order. A possible way of overcoming the memory and disk space limit, when 
finding the CCM equations, is to derive a perturbative solution to the CCM equations (as previously discussed in 
Sec. IV). We obtain the expansion coefficients of the ket-state correlation coefficients in terms of powers of some 
coupling parameter. As shown in Table III, we may then substitute the series expansions for the ket-state correlation 
coefficients into the ground-state energy equation, thus forming a power series for this quantity in powers of the 
coupling parameter. The prospective gain is that we might use less disk space in evaluating and storing the series 
coefficients than in storing the CCM equations. However, each CCM equation must be evaluated many times in order 
to determine the power series, which in turn means that we must significantly increase the speed of our code. We 
note that the number of LSUBm configurations that is needed in order to reproduce exact n th -order perturbation 
theory is a small fraction of the total number of LSUBm configurations. Hence, if one can state beforehand which 
configurations are necessary in order to reproduce n th -order perturbation theory, an elegant method of reducing the 
processing time and memory overheads might be to use this small fraction instead of the total number of LSUBm 
configurations. 

Another exciting prospect in order to increase the speed might be to generate the CCM equations in parallel. Also, 
the CCM is well suited to parallelisation as each CCM equation could be implemented on a separate processor. Note 
that the evaluation of each equation depends only on the target configuration and the form of the Hamiltonian, which 
remains constant. The possible weaknesses of this approach are that one is limited not only by the number of nodes 
in the system, but more importantly also by the "weakest link" . By this we mean that one can only progress as fast 
as the evaluation of the slowest equation at each iteration. However, very large gains are possible by sharing out the 
processing needs. 

To recap, we have mentioned that our results using the new CCM formalism are now fully competitive with QMC 
results, and we have seen that our approach is easily generalisable-|to, other systems. Possible futuee systems to which 
we may apply the new formalism include the valence-bond-solidsE2l'E3, systems with higher spincSEa, and models with 
electronic degrees of freedom such as the Hubbard modelEj. The nature of the ground states of new and interesting 
materials or spin models might be quickly and easily investigated by specifying the Hamiltonian, lattice and spin 
number, and development of the code might therefore lead to a powerful test-bench for various new ideas. The future 
also holds the possibility of very high-order calculations, which will increase our knowledge both of these systems and 
also of the CCM. With the inclusion of very high orders of approximation it is also hoped that the asymptotic nature 
of the CCM ground-state energies, sublattice magnetisations and phase transitions will become clear. In conclusion, 
we believe that the application of the CCM to the lattice spin systems is yielding excellent and interesting results, 
though there are still many avenues of fruitful and challenging research to be investigated. 
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FIG. 1. The square lattice LSUB4 bounding triangle is shown in this figure along with the LSUB4 lattice animals, illustrated 
by diagrams (a)-(e). The fundamental LSUB4 configurations for the planar model state are given by diagrams (a)-(j), and the 
fundamental configurations for the z-axis Neel model state form a subset of them, namely, all diagrams except (e), (i), and (j). 
The centres of the shaded squares mark the relative positions of the sites of the square lattice on which the spins are flipped 
with respect to the model state. 



17 




FIG. 2. Results for the CCM ground-state energy of the XXZ model on the 2D square lattice using the planar model state, 
compared to the Monte Carlo results of Ref. [44]. LSUBm critical A C1 and A C2 points are indicated by the boxes. 
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FIG. 3. Results for the CCM sublattice magnetisation of the XXZ model on the 2D square lattice using the planar model 
state. CCM results indicate non-zero, in-plane long-range order in the region — 1 < A < 1. 



19 



FIG. 4. The LSUB5 bounding triangle and symmetry-related correlation configurations for the triangular lattice. The centres 
of the shaded hexagons mark the relative position of the sites of the triangular lattice on which the spins are flipped with respect 
to the model state. See text for details. 
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(B) 
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5f vyt p tttttr 

(Mr torn to 

t f f r W Tr * 

FIG. 5. All 82 lattice animals of size 6 on a triangular lattice after symmetry reduction including translational and 12 
point-group symmetry operations (see text for details) . The centres of the hexagons mark the relative position of the sites of 
the triangular lattice on which the spins are nipped with respect to the model state. See text for a discussion of group A and 
group B diagrams. 
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FIG. 6. Results for the CCM ground-state energy for the triangular HAF. LSUBm critical points A C1 and A C2 are indicated 
by the boxes. 
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FIG. 8. Results for the CCM ground-state sublattice magnetisation for the triangular HAF. 
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TABLE I. Results obtained for the spin-i XXZ model on the 2D square lattice using CCM LSUBm approximations 
(m = 2, 4, 6, 8). Nf x denotes the number of fundamental configurations for the planar model state, which are further decom- 
posed in terms of connected and disconnected ones respectively, and Nf 2 denotes the number of fundamental configurations for 
the z-axis Neel model state. The ground-state energy per spin, E g /N, and the sublattice magnetisation, M + , at the isotropic 
Heisenberg point (A = 1) are shown, together with various critical anisotropy parameters. A C1 and A C2 indicate the LSUBm 
critical points for the planar model state corresponding to the ferromagnetic and antiferromagnetic phase transitions. A C3 
indicates the critical point for the z-axis Neel model state corresponding to the antiferromagnetic phase transition. Note that 
there are no terminating points in the LSUB2 approximation. 



m 


N Fl 


N F2 


EJN (A = 1) 


M+ (A = 1) 


A C1 


A C2 


A C3 


2 


1 (1+0) 


1 (1+0) 


-0.64833 


0.8414 








4 


10 (6+4) 


7 (5+2) 


-0.66366 


0.7648 


-1.249 


1.648 


0.577 


6 


131 (41+90) 


75 (29+46) 


-0.66700 


0.7273 


-1.083 


1.286 


0.7631 


8 


2793 (410+2383) 


1287 (259+1028) 


-0.66817 


0.7048 


? 


? 


0.8429 



TABLE II. Results obtained for the spin-i triangular-lattice HAF using CCM LSUBm approximations (m = 2, 3, 4, 5, 6, 7). 
Nf denotes the number of fundamental configurations which are further decomposed in terms of connected and disconnected 
ones respectively. Note that only CCM calculations up to the LSUB6 level of approximation are performed in this article. 
The ground-state energy per spin, E g /N, and the sublattice magnetisation, M + , at the isotropic Heisenberg point (A = 1) are 
shown for each LSUBm approximation. The terminating anisotropy parameters, A C1 and A C2 , which correspond respectively 
to a phase transition at A= — i and another believed [15] to be near A=l, are also given for each LSUBm approximation. Note 
that there is no terminating point in the LSUB2 approximation. 



m 


N F 


E g /N (A = 1) 


M+ (A = 1) 


A C1 


A C2 


2 


2 (2+0) 


-0.50290 


0.8578 






3 


8 (6+2) 


-0.51911 


0.8015 


-0.86 


5.47 


4 


30 (16+24) 


-0.53427 


0.7273 


-0.65 


2.20 


5 


143 (53+90) 


-0.53869 


0.6958 


-0.60 


1.98 


6 


758 (200+558) 


-0.54290 


0.6561 


-0.55 


1.77 


7 


4427 (837+3590) 


? 


? 


? 


? 
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TABLE III. Expansion coefficients in powers of A up to the 15th order for the ground-state energy per spin, E g /N, for the 
anisotropic spin-i triangular- lattice HAF obtained from the CCM equations in the LSUB6 approximation. The highest-order 
known exact series expansions up to the 11th order obtained by Singh and Huse [15] are also included for comparison. 
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U.UOll i uo 


n n^i ^4Q 

U.UO lOO^iiy 


8 


-0.0357291 


-0.0476598 


9 


0.0541263 


0.0685087 


10 


-0.0771681 


-0.1025446 


11 


0.1294578 


0.1565522 


12 


-0.1848858 
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13 


0.2857225 
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14 


-0.4463496 
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15 


0.7021061 
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